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Abstract. The analytical techniques of the Nekhoroshev theorem are used to 
provide estimates on the coefficient of Arnold diffusion along a particular resonance 
in the Hamiltonian model of Froeschle et al. (2000). A resonant normal form is 
constructed by a computer program and the size of its remainder ||i£ pt|| at the 
optimal order of normalization is calculated as a function of the small parameter e. 
We find that the diffusion coefficient scales as D oc \\Ro P t\\ 3 , while the size of the 
optimal remainder scales as ||-R op t|| <x ex P(l/ e °' 21 ) in the range f0" 4 < e < 10 -2 . 
A comparison is made with the numerical results of Lega et al. (2003) in the same 
model. 
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1. Introduction 

In a series of instructive papers, Froeschle et al. (2000, 2005), Lega 
et al. (2003), and Guzzo et al. (2005) presented the results of a de- 
tailed numerical investigation of the phenomenon of Arnold diffusion 
(Arnold 1964) in a Hamiltonian system of three degrees of freedom 
that satisfies sufficient conditions for the holding of the Nekhoroshev 
theorem (Nekhoroshev 1977, Benettin et al. 1985, Lochak 1992, Poshel 
1993). The aim of these studies was to establish quantitative estimates 
as regards a) the critical value of the small parameter e c below which 
the system enters into the so-called 'Nekhoroshev regime', and b) the 
dependence of the local diffusion coefficient D, along a particular res- 
onance, on e. In Guzzo et al. (2005) the authors reported to have 
also numerically observed global diffusion over an extended domain 
of the Arnold web. Other numerical studies over the years related 
to the same problem are: Kaneko and Konishi (1989), Wood et al. 
(1990), Dumas and Laskar (1993), Laskar (1993), Skokos et al. (1997), 
Giordano and Cincotta (2004). Some early results of the literature are 
briefly commented in the discussion. 

In what follows we focus on an analysis of point (b) above. According 
to Lega et al. (2003), a tedious numerical calculation yields that the 
local diffusion coefficient D along a resonance depends monotonically 
on e. While the limited range of the numerical data did not allow for 
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a precise fitting, the authors presented evidence that the numerical 
function D(e) was decreasing faster than a power-law D(e) oc e b , while 
it was consistent with an exponential law D{e) oc exp(— l/e a ) for some 
unspecified constant a. A subsequent study in a mapping model anal- 
ogous to the above Hamiltonian model yielded the concrete estimate 
a ~ 0.28 (Froeschle et al. 2005). 

In the present paper our aim is to provide a deeper understanding of 
these numerical results by having recourse, precisely, to the analytical 
apparatus of the Nekhoroshev theory, that is, the construction of a 
resonant Birkhoff normal form and of estimates on the size of the 
remainder R of the normal form at the optimal order of normaliza- 
tion. The optimal remainder function R opt is a crucial quantity for the 
dynamics because it represents the true size of the perturbation of the 
system that corresponds to its deviation from an integrable system. 
Indeed, in Nekhoroshev theory all bounds on the variations of the 
actions follow from estimates on the size of R op t- In order, therefore, to 
establish the link between the Nekhoroshev theory, on the one hand, 
and the problem of the speed of diffusion, on the other hand, in the 
present paper we seek to determine: 

a) the dependence of the local diffusion coefficient D on the size of 
the remainder ||.R op i|| at the optimal order of normalization, and 

b) the precise dependence of ||-R pt|| on e. In particular, we seek to 
calculate ||i? pi|| when e becomes small enough so that the exponen- 
tially small scaling of ||-R pt|| with e clearly shows up. This required 
to consider values of e one order of magnitude smaller than the values 
considered in the numerical experiments of Lega et al. (2003). 

We can immediately state a brief summary of our main results. We 
find that: 

a) D scales with ||-R pt|l as a power law D oc ||i? opt || c , with c ~ 3. 
To determine this relation, use was made of the data of Lega et al. 
(2003) for the numerical values of D, and a program was written by 
the author in order to compute the Birkhoff normal form and thereby 
the remainder R op t- 

b) In the range 10~ 4 < e < 10~ 2 , ||-R op t|| scales as ||-R op t|| = 
exp(— b/e a ) with a ~ 0.21 and b ~ 3.22. This also yields the exponential 
scaling D oc exp(— 36/e a ). Note that the exponentially small scaling 
shows up clearly only for e < 0.001, while in the range 0.001 < e < 0.02 
a power-law fitting 1 1 R op t \ \ °^ e 2 ' 45 yields marginally better results than 
the exponential fitting. In fact, we find that although the onset of the 
'Nekhoroshev regime' can be placed at a threshold value e ~ 0.01, the 
deviation of the remainder from a power law becomes clear only for e 
one order of magnitude smaller than this value. 
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The computational algorithm used to obtain the above results is 
described in section 2. It follows essentially the normalization scheme 
described in the lecture notes on exponential stability of Giorgilli (2002, 
pp. 86-87). This scheme is different from the traditional scheme based on 
normalization by successive orders of e. We further modify it to avoid 
using K— truncated resonant modules (see section 2). People working 
on mathematical aspects of the Nekhoroshev theory are definitely famil- 
iar with all related notions. However, a somewhat extended description 
is given in section 2, to help rendering the subject more accessible also 
to people oriented towards the applications. 

Section 3 then focuses on the implementation of the algorithm in the 
particular model of Froeschle et al. (2000). By constructing the normal 
form and identifying the optimal order of normalization, the size of 
R op t is evaluated for various values of e. This yields the exponential 
scaling of ||-R op t|| on e. Using also the data of Lega et al. (2003) yields 
the power-law relation of D with ||i? opt ||. Section 4 contains a brief 
discussion of the results. 



2. The normal form algorithm 

2.1. The classical construction 

The analytical part of the Nekhoroshev theorem is based on the con- 
struction of a normal form for a n degrees of freedom Hamiltonian 
system of the form 

H(J,0)=H o (J)+eH 1 (J, ( f > ) (1) 

where ( J, (p) = (Ji, J2, J n , <f>i, <f>2, ■■■4>n) are action-angle variables, 
Hq satisfies appropriate non-degeneracy and convexity conditions, and 
H is analytic in a complexified domain of the actions and the angles. 
We shall be concerned with a local construction of a normal form for the 
Hamiltonian (1), valid within some (small) preselected open domain W a 
of the action space (the index a means action space). The analyticity 
condition implies that if Hi is Fourier expanded 

H 1 (J,<f>)= £ H ltk ( J) exp(ik ■</>)) (2) 

there are positive constants A, a such that for all values of J G W a the 
coefficients Hij, are bounded by 

\H hk (J)\ < Aexp(-\k\a) . (3) 
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The function Hi is analytic in the domain W a x T a n , where T a n = 
{(f) : Re0 £ T™, |Im</>j| < o~,i = l,...n}. Furthermore, bounds on the 
function iii can be found in domains W a x TV™, for any positive <r' < cr. 

Now, the purpose of the normal form construction is to perform a se- 
ries of canonical transformations (J, <f>) — ► (J*- 1 -*, 0^), — ► (</ < - 2 \ f^ 2 -*); • • •> 
such that, after r normalization steps, the Hamiltonian is expressed as 

#M(j(r)^(r)) = H{J{j(' r \^)<P{j( r \^)) 

= Z (r) (J (r) ,</> (r) ) +i? (r) (J (r) ,0 (r) ) . (4) 

In the last expression, the first term Z( r \j( r \ (f)^), called the normal 
form, is constructed in such a way as to yield a simple dynamics in the 
variables (J^ r \(f)^) (e.g. in the non-resonant or simply resonant case 
the Hamiltonian Z( r \j( r \ <f>^) is integrable). On the other hand, the 
second term R( r \j( r \ 4>^), called the remainder, represents a pertur- 
bation to the dynamics induced by the normal form. The goal of the 
theory is to proceed with the normalization as far as necessary in order 
to minimize the size of the remainder. In exceptional cases of integrable 
Hamiltonian systems one may proceed by infinitely many steps, within 
some domain of convergence, and show that the remainder tends to 
zero in the limit r — > oo. In the generic case of non- integrable systems, 
however, one can only reduce R( r \j( r \ (f)^) to a finite minimum size. 
This is found at a specific order of normalization which is hereafter 
called the optimal order of normalization r op t . The size of the remainder 
i?( r )(j( r ) 5 </>( r )) for any other order r ^ r op t is bigger than the size of 
i?( r ° p *)(j( r ), (f>^). Standard theory yields estimates r opt ~ l/e b and 
\\R(ro P t)\\ „ exp(-l/e a ) for some positive integers a, b depending on 
the number of degrees of freedom. 

The classical construction leading to expressions of the form (4) is 
nowadays analyzed in many books (see e.g. Boccaletti and Pucacco 
(1996), Contopoulos (2002), Morbidelli (2002), or Ferraz-Mello (2007) 
for instructive introductions). The main steps can be summarized as 
follows: 

a) One makes a choice as regards which Fourier terms Hi^(J) exp(ifc- 
0) in Hi will be retained and which terms will be eliminated in the 
normalized Hamiltonian. The terms to be retained are specified by their 
integer vector k belonging to a particular subset of 7L n called resonant 
module (denoted hereafter by M)). The choice of resonant module is 
dictated by the location of the domain W a in the action space, i.e., by 
whether this domain is close to or far from particular resonances (see 
subsection 2.3). 

b) The Fourier terms not belonging to M are the ones to be elim- 
inated by a canonical transformation. Let us consider the first nor- 
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malization step. We shall consider canonical transformations of old 
variables (J, (f>) to new variables (J*- 1 ^, (f>^) obtained via a Lie generat- 
ing function xi(J^\ <fi^)- The function Xi(J ( - 1 ^, 4>^) will be selected 
so as to be of the 'first order of smallness' (the precise meaning of this 
sentence is analyzed below). The canonical transformation is defined 
by 

J = exp(L Xl )J (1 \ 4> = eML Xl )<t> {1) (5) 

where L X1 = {-,xi} is the Poisson bracket operator, and exp(L Xl ) is 
the formal exponential development of L X1 . A well known property of 
Lie transformations is that for any function of the canonical variables 
f(q,p) one has 

exp(L x )f(q,p) = /(exp(L x )g,exp(L x )p) . (6) 

As a result, the transformed Hamiltonian, after the action of the gen- 
erating function \i is given by 

= exp(L Xl )H = H + L X1 H + ±L 2 X1 H + . . . (7) 

In order to specify the function Xi{J^\ 4^)i we select from H\ all 
the terms satisfying the following two conditions: a) being of the 'first 
order of smallness', and b) not belonging to M. Denoting by h\ the 
sum of these terms, the terms are eliminated by selecting xi so as to 
satisfy: 

L Xl H + hi = . (8) 

Equation (8) is called homological. It is reduced to a diagonal system 
of linear algebraic equations. Namely, writing hi as 

hi = hi :k (J) exp(ifc • 4>) 
k 

and setting 

xi = ^2xi,k( J ) ex P( ik • 4>) 
k 

we find, through (8), the solution for the coefficients Xi,k(J) given by: 

/ n hik(J) 

X ^ J) = ikMT) (9) 

where 

cj(J) = VjHo(J) (10) 

is the frequency vector of the unperturbed Hamiltonian at the point J 
of the action space. 
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c) In the same way we eliminate the terms of the second order of 
smallness in or, in general, the terms of the r-th order of smallness 
in 

jy(r-i). The 

recursive formulae for the r-th canonical transformation 
J(r-i) = eX p(L Xr )j( r ), 0( r " 1 ) = exp(L Xr )0( r ) are: 

{H ,Xr} + 4 r ~V =0 

fl-W = exp(L Xr )#( r - 1 > . (11) 

Following iteratively the above procedure, the function takes the 
form (4). The remainder consists of all the terms of of 'order of 
smallness' larger than r. It can be written as 

flW = h% + h% + ... (12) 

This series is an analytic function in a restriction of the domain of 
analyticity of the original hamiltonian provided that the Lie transfor- 
mation in (11) is convergent at every step (which is ensured for 'small 
enough' generating functions Xr)- 

2.2. Book-keeping 

A clarification of the meaning of 'order of smallness' is now in order. 
When one implements the recursive scheme of Eqs.(5) and (11), one 
needs to decide how should at every step be split into terms of 
different orders. Such a splitting should be based on the size of each 
term relative to the size of the other terms. In a computer program, it 
is customary to introduce a so-called 'book-keeping' parameter A. This 
is a parameter the powers of which appear as coefficients in front of 
the various terms in the expansion. A term with coefficient X p is said 
to be 'of order of smallness' p. This helps separating the terms in the 
program by, say, recalling the Coefficient(/i (r_1 ), A, r) in an algebraic 

program like Mathematica in order to get the term hr' ^ in Eq.(ll). As 
the program carries on the normalization, the book-keeping coefficients 
change according to the rule that the coefficient of the Poisson bracket 
of two terms with coefficients X p and X q is X p+q . In the end we set A = 1 
to recover the correct values of all the coefficients. 

This way of organizing the series is practical, but it also reflects 
a fundamental decision on what quantity we consider to represent 
'smallness' in the perturbation series. If the overall size of H\ is an 0(1) 
quantity, then it is customary to identify e itself as the small parameter. 
In that case the normalization scheme proceeds by ascending powers of 
e, and the generating functions Xr are of order 0(X r ) = 0(e r ). We shall 
see in subsection 2.4 how to modify the 'book-keeping' so as to separate 
terms of different smallness which co-exist within H\ , and then within 
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H^ r \ provided that the function Hi satisfies analyticity conditions of 
the form (3). 

2.3. Resonant K-module 

The choice of resonant module depends on the location of the domain 
W a in the action space and on the properties of the unperturbed Hamil- 
tonian Ho. Let us consider the homological equation of the first order 
under the choice of book-keeping A = e. This reads: 

{H , X i} = -eH 1 (13) 

where H\ is the part of H\ containing terms to be eliminated. The 
solution of (13) is 

H hk (J)exp(ik ■ <j>) 

Xl = \ h» ik^J(J) • (M) 

Provided that the resonant manifolds k ■ uj( J) = are dense in the 
action space, for any small open domain W a there are infinitely many 
resonant manifolds crossing W a - This implies that there is a dense set 
of values of the actions J within W a for which a divisor in (14) becomes 
exactly equal to zero. In practice, this means that a generating function 
like xi i n (14) cannot in fact be determined in any open domain W a . In 
Nekhoroshev theory, an often stated remedy to this problem consists of 
splitting the set of Fourier harmonics exp(ik-(f>) into low order and high 
order harmonics according to whether \k\ = |fci| + . . . + \k n \ is smaller 
or larger than some positive integer K (see, for example, Morbidelli 
and Guzzo 1997). The choice of value of K is arbitrary, but the choice 
K ~ l/e a , for some constant a, is suggested by the request that there 
be no overlapping of the resonant domains (of width 0{e 1 / 2 )) formed 
around each of the resonant manifolds k-co(J) = 0, with \k\ < K, other 
than the overlapping near the loci at which the manifolds intersect 
each other (see Morbidelli 2002, p. 119). Fixing K and requesting that 
only the harmonics with \k\ < K be normalized ensures that divisors 
k ■ uj(J) with \k\ > K will not appear in the series. Furthermore, the 
domain W a is crossed by only a finite number of resonant manifolds 
with |fe| < K. More precisely, denoting by lZ a the subset of integer 
vectors k for which the resonant manifolds k ■ u(J) = intersect the 
domain W a , and by bn a a subset of lZ a such that any vector of 7Z a 
can be produced by a linear combination of the vectors of bn a with 
integer coefficients, the relevant property is that bn a contains only a 
finite number of vectors k with \k\ < K and an infinite number of 
vectors with \k\ > K. Denoting by b^ a the subset of vectors k G bn a 
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with \k\ < K, the resonant module M. a ,<K is determined by: 

M a ,<K = span{6|f } . (15) 

Let M a ,>K be the complement of A4 a ,<K with respect to lZ a . Under 
the 'book-keeping' choice A = e, one then has the following iterative 
algorithm of determination of the normal form locally, i.e., within the 
domain W a : 

a) Assuming that (r — 1) normalization steps were completed, split 
Hr ^ as Hr ^ + h!f 1 \ where Hr" ^ contains all the terms with 
\k\ < K not belonging to M a ,<K- 

b) Determine the generating function Xr by solving 

{H ,xr} = -H^ l) (16) 

c) Find the new Hamiltonian 

= exp^J^-i) = Z +eZ 1 + . . .+e r Z r +H^ s +H$ nres (17) 

in which Zj, j = 0, . . . , r are normal form terms belonging to the reso- 

(r) 

nant module M a ,<K, Hr es contains all the Fourier terms belonging to 

(r) 

•Ma,>K, and Unonres contains the remaining Fourier terms that are not 
in any of the previously determined sets. This completes one step of the 
iteration algorithm. As shown in Morbidelli and Giorgilli (1997), while 
all the terms of order higher than r in these two functions are bounded 
by a quantity 0(exp(— Ka)), i.e., they are both exponentially small in 

(r) 

K, the terms of Hres are the ones making the leading contribution to 
the remainder. 

2.4. A MODIFIED ALGORITHM 

Two problems arise in the implementation of the above algorithm: 

a) Eqs.(9) or (14) imply that the denominators of all the series 
coefficients are functions of the actions, in the form of products of 
divisors k ■ lo(J). This implies that even if Hq{J) is polynomial in the 
actions, the Fourier coefficients of the terms exp(ifc • 4>) in or %( r ) 
are in general rational in the actions. The extra computational load of 
dealing with rational expressions (and their derivatives entering into 
Poisson brackets) makes the whole algorithm hardly tractable in this 
form. We defer the solution of this problem after the analysis of point 
(b) below, to which it is linked. 

b) One needs to store many more intermediate coefficients than those 
eventually needed in order to have a complete knowledge of the trans- 
formed Hamiltonian or of the generating function within any specified 
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domain of interest. This is a rather technical problem, the solution of 
which is, nevertheless, crucial in constructing an efficient algorithm for 
the computation of the normal form. A detailed quantitative discussion 
of this problem is deferred to the appendix. Here we state the main 
result: Under the book-keeping scheme A = e, assume we are interested 
in specifying the transformed Hamiltonian in a 'domain of interest' 
defined by the maximum orders (r max , K max ) , where r max denotes the 
maximum order in the book-keeping variable e and K max denotes the 
maximum order in Fourier space. Then, in order to be able to specify 
all the terms belonging to the 'domain of interest', we must store, at 
intermediate normalization orders r < r max , all the Hamiltonian terms 
of Fourier order \k\ < K(r max — r) + K max . As shown in the appendix 
these intermediate terms outnumber by far the terms finally stored 
within the 'domain of interest'. 

Issue (b) can be resolved by introducing a different book-keeping of 
the terms, on the basis of their Fourier order rather than e— order. This 
organization of the series is analyzed in detail in Giorgilli (2002, pp. 86- 
87) and we present the main points of this analysis in the appendix. 
Essentially, it reflects the fact that for any order r < r max , many Fourier 
terms with a coefficient e r in front, which would thus be formally stored 
as of order A r under the book-keeping e = A, are actually of much 
smaller size than e r , because the size of any exp(i/c • <j>) term scales as 
(e — CT )l fe l by virtue of the analyticity condition. This scaling is already 
present in the original Hamiltonian and it propagates at all subsequent 
normalization steps. Precisely, this fact is recognized by 'book-keeping' 
the terms on the basis of their Fourier rather than e-order. The exact 
algorithm, which replaces the algorithm of subsection (2.3), reads as 
follows: 

1) Define K' = max{[l/cr], 1}. 

2) Place a book-keeping coefficient A p in front of each Fourier term 
of the form exp(ifc • (f>) in Hi, where p = [\k\/K'] + 1. 

3) Carry on the normalization (Eqs.(ll)) by successive powers of A. 
As shown in the appendix, with such an algorithm there are no extra 

terms, outside the domain of interest, that need to be computed. 

Furthermore, it is possible to incorporate a solution to issue (a) in 
the same algorithm. First, we notice that, as defined in step (1) above, 
the constant K' does not pose an upper bound (like K) in the order of 
the Fourier harmonics being normalized at successive steps, since the 
algorithm requests that terms of increasing order rK' be normalized at 
the r-th normalization step. In principle, this would cause a problem 
when r becomes larger than K/K' . But in practice, one wishes to avoid 
the problem of the appearance of the actions in the denominators of 
the Fourier coefficients for both \k\ < K or \k\ > K, that is for both 
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r < K/K' and r > K/K' . It can be shown (section 3) that if the 
domain W a is 'resonant', that is M a: <K / {0, 0, ... , 0}, the width of 
W a scales with e as 0(e 1 ^ 2 ). Non-resonant domains can also be covered 
by subdomains of size 0(e l l 2 ). Choosing a point J* in the interior of 
the domain W Q , we then develop the Hamiltonian locally, within W a , 
by expanding the actions as 

J = J ++e 1 /2 / _ ( 18 ) 

The transformation J — > I is not canonical under the usual Poisson 
structure, but the dynamics remains unaltered if a new Hamiltonian 

H' = e- 1/2 H (19) 

substitutes H as generator of the equations of motion. For example, in 
a Thirring-type model 

j2 , j2 

H = 1 - 2 + e Y, c kiM exp(i(fci^i + fe^)) (20) 

fel,fe2 

we get, introducing also the 'book-keeping' parameter A: 
i?' = const + Ji*/i + J 2 ,/ 2 + g 1/2 fA x ^ 2 

+ ^ A[l fc l^'] +1 c fcl , fc2 exp(i(A; 1( Ai + A; 2( A 2 ))) . (21) 

Consequently, the perturbation has been rescaled to e 1//2 , i.e., it follows 
the size 0(e 1 ^ 2 ) of the domain W a . Nevertheless, the unperturbed 
frequencies are now lu(J*) = thus they do not depend on the 
actions, since the terms quadratic in the actions are now formally 
of order A. This implies that the divisors are of the form k ■ lj*, i.e., 
independent of the actions. Consequently, the algorithm proceeds with 
polynomial rather than rational (in the actions) Fourier coefficients at 
every step. This also means that there is no longer need for introducing 
a K-truncated resonant module. Thus, in all calculations the resonant 
module is specified only by the resonances between the frequencies 
uj*. If lj* are incommensurable, one may specify a module correspond- 
ing either to no resonance or to an approximate resonance between 
the frequencies w*. As discussed by Morbidelli (2002, p. 48-49), such a 
construction is still consistent with the appearance of an exponentially 
small remainder at the optimal order of normalization. This is because, 
as demonstrated in the appendix, the order of normalization has a 
linear relation with the maximum Fourier order of the terms contained 
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in the normal form. Furthermore, we find that r opt ~ e -1 / 2 which also 
implies \k\ op t ~ e -1 / 2 . Hence the dominant terms in the remainder are 
of Fourier order |fe| > \k\ op t, i.e., they have a size 0(exp(— ae~ 1 / 2 )). 
This means that \k\ op t plays now a role similar to K in subsection 2.3. 
It is interesting to note, however, that such a result is obtained here 
a posteriori, i.e., after the construction of the series. That is, with the 
present algorithm one may proceed by calculating the series up to any 
desired value of r, and then by checking whether the optimal order 
r opt was reached. This will also specify automatically the optimal order 
|fe|opt- On the contrary, following the algorithm of subsection (2.3), one 
would have to fix a value of K in advance and then check whether such 
a value yields the optimal estimate for the remainder. Otherwise, the 
calculation should be repeated from the start, by making a different 
choice of K. 



3. Results 

3.1. Hamiltonian model and resonant dynamics 
The Hamiltonian model of Froeschle et al. (2000) reads: 

H = H + eHi = 1 2 + h + (22) 

2 4 + COS (f>l + COS <p2 + COS (f>3 

where (Ii,cj)i), i = 1, 2, 3 are action - angle variables and e is a pertur- 
bation parameter (note the change of notation for the actions, J I, 
with respect to the previous section, for consistency with the notation 
of Froeschle et al. (2000)). 

The integrable part Hq produces a simple dynamics: Ii = and 
4>i = ojo,i = h, 4>2 = ^0,2 = h, <j>3 = ^o,3 = 1- Thus all three actions 
are integrals of motion and all three angles circulate with constant 
frequencies. The surface of constant energy H = E is a paraboloid in 
the action space given by if + if + 2/3 = 2E (Fig. la). On the other 
hand, the resonances fcia>o,i + ^2^0,2 + ^3^0,3 = k\I\ + kill + k^ = 0, 
with k = (^1,^2,^3) £ Z 3 , \k\ = I fei| + |fe2| + |fes| / 0, define planes 
in the action space which are always normal to the plane 12)- The 
intersections of the resonant planes with the surface of constant energy 
define the web of resonances on this surface which is a set of parabolic 
curves (Fig. la). When viewed from the top of Fig. la, the projection of 
these curves on the plane (Ji, I2) defines a set of straight lines (Fig. lb). 

When e = all the points on the surface of constant energy of Fig. la 
are neutral equilibria, which correspond to invariant 3-tori in the six- 
dimensional phase space. If, however, e / 0, an O(e) volume of invariant 
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Figure 1. (a) The white grid-lined surface is a part of the paraboloid of constant 
energy l\ + 7| + 2/3 = 2E in the action space for the value of the energy E = 1. 
The three gray-shaded planes correspond to the resonances Ji — 2/2 = (parallel 
to the double arrow), 3Ji + I2 — 1 = 0, (to the left of the double arrow), and 
2/i + 2/2 — 1 = (to the right of the double arrow). The intersection of all the 
resonant planes with the paraboloid of constant energy produces a set of parabolic 
curves which is the 'web of resonances', (b) The projection of the web of resonances 
kill + /C2/2 + &3 = for \k\ < 5 on the (Ji, J2) plane. Whenever the coefficient of 
the term (for one particular resonant vector k) of the Fourier development of Hi is 
positive, we only plot a single line corresponding to the projection of the associated 
parabolic curve on the plane (7i, 12). If this coefficient is negative, we plot the central 
line and two other parallel lines marking the borders of the resonance as obtained by 
an analytic estimate of the associated separatrix width when e = 0.003. These rules 
correspond to a Poincare surface of section <j)R = 0, where (f>n is the resonant angle 
associated with each resonance (see text for details), and due to this reason Fig. lb 
is directly comparable to Fig.l of Lega et al. (2003). The analysis in that paper and 
below refers to orbits exhibiting chaotic diffusion along the direction indicated by 
the double arrows at the borders of the resonance Ii — 2h = 0. 
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tori of the phase space are destroyed, according to the KAM theorem, 
and replaced by a chaotic subset of orbits. The constant energy condi- 
tion foliates the phase space into 5-dimensional hypersurfaces. For any 
fixed value of the angles, the projection of one hypersurface on the space 
of actions defines a manifold resembling to a paraboloid, like in Fig. la, 
which, however, has some deformation of order 0(e). The chaotic orbits 
can drift on this manifold. We shall focus on orbits sliding on the 
manifold of constant energy along the chaotic borders of resonances 
in the direction indicated by the arrows of Figs.la,b. As shown below, 
the planes of resonances k ■ I2, 1) = of the unperturbed system 
are transformed to resonant manifolds of the perturbed system that 
also resemble to planes, but with some thickness of order 0(\hke\ 1 ^ 2 )), 
where hk is the coefficient of the term exp(ik ■ <j>) in the Fourier devel- 
opment of the perturbation H\. Thus the intersections of the resonant 
manifolds with the manifold of constant energy form resonant zones 
of thickness 0(|/jfce| 1//2 ). When projected on the plane (h,h) the web 
of resonances looks like in Fig. lb (the reason why in some cases we 
show one single line for one resonance while in other cases we show a 
triple line is discussed below) . The chaotic orbits studied by Lega et al. 
(2003) are calculated at a value of e at which the system is in the so- 
called 'Nekhoroshev regime', i.e. the width of the resonant zones is small 
and there is no significant resonance overlap. The orbits drift slowly 
along the border of the resonance (k\, tt2, ks) = (1, —2, 0), starting with 
initial conditions close to, but outside the intersection of this resonance 
with a different resonance, namely (k±, k2, ^3) = (3, 1, —1). The chaotic 
drift along the zone of the resonance (1,-2,0) produces a slow secular 
change of the value of the action Ip = 2ii + I2 (see Eqs.(26 - 28) 
below), which is the action associated with the oscillation in the so- 
called "direction of fast drift", i.e., normal to the resonant plane. The 
change of the value of Ip with time was used by Lega et al. (2003) in 
order to measure numerically the coefficient of the chaotic diffusion. 
The latter causes also a slow change of the value of the action ^3, due 
to the confinement of the orbits on the manifold of constant energy. 

An estimate of the width of resonances can be made as follows. If 
we complexify the angles in a domain of C containing the real 3-torus, 
i.e., if we make the replacement 4>j — ► (f)j + iuj in Hi, with <f>j,Uj real, 
then in any direction of the space (u±, 1x2,1x3) the quantities Uj can 
vary from zero up to a value at which Hi becomes singular. However, 
there is a lower bound on the distance of all the singularities of Hi 
from the real 3-torus, i.e., from the value (ui, 112,113) = (0,0,0). Thus, 
the values of Uj at the point of the closest singularity cannot be all 
infinitesimally small. In fact, if we develop the denominator of Hi up 
to terms of second degree in Uj, we find an estimate of the position 
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Figure 2. The moduli of the coefficients hk of the Fourier development (23) as a 
function of the modulus \k\ for all resonant terms up to \k\ — 40. The straight solid 
line is an exponential fit to the maxima of the \hk\ values as a function of |fc|. 




of the singularities of Hi close to the real 3-torus given by the values 
Uj = Uj jS satisfying 



u 



l s cos 4>i + u\ s COS 02 + u% s COS 03 + 8 + 2(cos 01 + COS 02 + COS 03) = 
and 

iti )S sin 0i + u 2 , s sin 2 + u 3jS sin 3 = . 
The first of the above equalities implies 



1,5 + «2,a + U 3,s > Wi, s COS 4>l + u\ s COS 2 + u\ s COS 3 | > 2 



u 

The interior of any parallelepiped in the u-space defined by the three 
inequalities \uj\ < Gj, for some positive constant numbers (Tj, j = 1, 2, 3, 
constitutes a domain of analyticity of Hi if it is contained within the 
sphere u\ + u% + = y/2. The optimal estimate for the domain of 
analyticity corresponds to the parallelepiped with maximum volume 
8<7i<T2CT3, that is to the choice o\ = a 2 = C3 = \/2/3 = 0.82.. ~ a. 

On the other hand, a more precise value of a can be found by 
expressing Hi triple Fourier series 

-. oo oo oo 

= y y y h k eMik ■ 0) (23) 

4 + COS 01 + COS 02 + COS 03 

where k = (ki, k 2 , k%), = (0i, 02, 03). According to the Fourier theo- 
rem on analytic functions, the coefficients \hk\ are bounded from above 
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by 

\h k \ < Aexp(-|A;|cr) (24) 

with \k\ = \k\\ + \k2\ + l^l and A, a positive constants. Thus Eq.(24) 
predicts that the coefficients h k decay exponentially. This is shown in 
Fig. 2, where we plot the values \h k \, calculated by a series development 
of H\ using Mathematica up to the 40th order, against \k\. A numerical 
fitting to the upper bound of this diagram (solid line) yields A = 0.046, 
and a = 0.87. The latter value is close to the estimate a = a/2/3 = 0.82 
found above heuristically. 

Having determined the size of the Fourier coefficients h k , the normal 
form construction described in subsection (2.4) eliminates all harmonics 
in the Hamiltonian perturbation H\ (written in the form (23)), except 
for the harmonic h k exp(ik-(fi) corresponding to the particular resonance 
into consideration. To the lowest order, the resonant normal form reads: 

p , j2 

H res = 1 2 + h + 2e% cos(fc • <P) + ... , (25) 

where we have used the property h k = h- k following from the reality 
and even parity with respect to the angles fa of the Hamiltonian per- 
turbation H\ . Introducing resonant action - angle variables (I R ,(j) R ), 
(If, 4>f) via the generating function: 

S = (k ■ 4>){I R - k l ) + (mi^i + m 2 fc)I F + fol's (26) 

where the integers m\,m2 are any pair satisfying m\k\ -\-m2k2 = 0, the 
resulting canonical transformation is 



^3 ^3 

kl(lR ~ 2 2 ) + m^F, I 2 = k 2 (lR ~ 2 2 ) + m 2^F, 

<Pr = hfa + k 2 (p2 + k 3 <p3, <t>F = mi^i + m 2 (f>2, (27) 

^3 = I3, 03 = <f>3 ■ 

By virtue of (27), the resonant Hamiltonian (25) takes the form (apart 
from a constant): 

H res = \{m\ + m 2 2 )I F + I 3 + ^(kf + kl)I 2 R + 2eh k cos <f> R (28) 

i.e., it is split in two parts depending only on the actions If, h, and a 
third part which is the pendulum Hamiltonian: 

H pend = ^{kl + kl)I 2 R + 2eh k cos<t ) R . (29) 
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The actions ^3 and Ip = (rn\I\ + 771-2/2) / '{m\ + rn^) are integrals of the 
Hamiltonian H res (in the case of the resonance (k±, k2, k%) = (1, —2, 0) 
we have m\ = 2, 777,2 = 1 and the precise definition of the fast action is 
If = (2/i +/2)/5 which differs by a factor 1/5 from the definition given 
at the beginning of this section for that particular resonance). On the 
other hand, the width of the resonance is determined by the separatrix 
half-width of the pendulum Hamiltonian H penc i 



The phase portrait of H pen d is shown schematically in Fig. 3. In reality, 
the ideal separatrix given by (29) should be replaced by a thin chaotic 
layer produced by the weakly chaotic motion near the separatrix due 
to to higher order coupling terms of the Hamiltonian. When e is small, 
however, the width of the chaotic zones is very small. In that case, 
Eq.(30) can be used to approximate the maximum normal distance to 
the resonance between the upper and lower branches of the separatrix- 
like chaotic zones. 

The resonant action Ir changes in time according to 



Any variation AIr of the resonant action results in variations AIi, 
A/2, such as to respect the integral Ip, i.e., Alp = 0. Eq.(27) then 
gives 



When the resonance web is projected on the (Ji,/2) plane, as in 
Fig. lb, the variations AIi and A/2 of Eq.(31) correspond to motions 
across the resonance, i.e., along a line normal to the resonance line 
k\Ii + k2h + ^3 = (this is called the "direction of the fast drift" by 
Froeschle et al. (2000)) On the other hand, when the web of resonances 
is visualized numerically in the action space, as, for example, by the 
method of the Fast Lyapunov indicator (Froeschle et al. 2000), one 
has to make an appropriate choice of Poincare surface of section in 
order to eliminate the angles from the calculation. The choice ^3 = 0, 
\<f>i I + 1^2 1 < 0.05 made by these authors corresponds essentially to 
setting the resonant angle <pR = k\4>i + k2<t>2 + &303 to a value very 
close to zero, i.e., <pp ~ 0. We then see that this means to plot (a) two 
sets of points in the plane (ii,^), corresponding to passing close to 
the maxima or minima of the theoretical separatrices, when < 0, 
or (b) one set of points on the plane (/i,^), corresponding to passing 
close to the X-point of each separatrix, when hk > 0. Consequently, 




(30) 



I R = -dH pend /d(f)R = eh k sm(fi R 



Ah = kiAI R , AI 2 = k 2 AI R 



(31) 
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Figure 3. The pendulum phase portrait for H pen d given by Eq.(29) (arbitrary units). 
The left and right dashed vertical lines indicate the section 4>r = when hk > 
and hk < respectively (i.e. <pr = (f>R and ifiR — (f>R — ir respectively). 

the thin chaotic borders of the resonances will appear as follows on 
the surface of section: In case (a) (hk < 0) we see two lines parallel 
to the resonant line ^1/1 + ^2/2 + ^3 = 0. These lines define a zone of 
width determined by Eq.(31). In case (b) (hk > 0) we see only one line 
coinciding with the resonance line itself. These rules are followed in the 
plot of Fig. lb, for all the resonances with \k\ < 5. This figure is to be 
compared with Fig.l of Lega et al. (2003). The two figures compare 
well not only qualitatively, but also quantitatively, i.e,, the theoretical 
resonance widths found above are very close to the widths determined 
numerically by the FLI method. 

3.2. Normal form and Nekhoroshev estimates 

The resonance dealt with by Lega et al. (2003) is I\ — 2I2 = 0. The dif- 
fusion takes place on the thin chaotic border along this resonance (Fig. 2 
of Lega et al. (2003)) when the initial conditions are taken in a small 
region of the thin chaotic border. In order to perform a transformation 
of the type given by Eqs.(18) and (19), we have to specify the central 
values (ii*,/2*) = (^i*,W2*) with respect to which the Hamiltonian is 
expanded. By visual inspection of Fig. 2 of Lega et al.(2003) a conve- 
nient choice compromising the central values of in all the panels is 
Iu = 0.31, I2* = 0.155. Renaming the variables of Eq.(22) by the same 
symbols according to e~ 1//2 i? — ► H, e^ x l 2 (I{ — 1^) — > ij, i = 1, 2, 3, and 
Fourier-expanding up to order 40, the Hamiltonian reads (apart from 
a constant): 

H = H *(I) + e 1/2 H u (I, <j>) = uj u h + UJ2J2 + h + 
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eV*( I -l±l+ £ c k eMik-®) . (32) 

V |fc|<40 ' 

According to the definition of K' given in subsection 2.4, the term 
e l l 2 (ll + /|)/2 should be considered of order of smallness p = 1. Then, 
it is convenient to slightly modify the definition of K' so as to render 
e~ K a comparable to the factor e 1 / 2 . We thus set: 

K , r iog(^ 1/2 ) 

The range of values of e for which Lega et al. (2003) provide numerical 
data on the diffusion coefficient is 30 < e _1 < 1000. For most values 
within this range K' is equal to a constant value K' = 3. We thus 
simply set K' = 3 in all the calculations. This means that the Fourier 
modes of H\* with < \k\ < 2 are book-kept as of order 1, modes 
with 3 < \k\ < 5 as of order 2, etc. The terms of Ho* are of order zero, 
i.e., they are the ones to be used in the kernel operator {-,ifo*} of the 
homological equation. 

Although the computer program allows the user to carry along the 
powers of e 1 / 2 in the expansion, the extra requirement of memory in 
order to store the associated exponents of e 1 / 2 is prohibitive. We thus 
had to make separate runs of the computation of the normal form 
for various numerical values of e. That is, one numerical value of e 
was supplied to the file storing the coefficients of (32) for each run. 
Even so, a memory limit of 2GB was reached by computing about 
5 x 10 7 complex coefficients of the Hamiltonian and of the Lie generating 
function per run. Such a number of terms corresponds to a truncation 
r < 15 in the normalization order, or K = 3 • 15 — 1 = 44 in Fourier 
space. Despite these limits, it was in the end possible to observe the 
asymptotic properties of the Birkhoff series within the whole range of 
values of e down to a value e = 0.0001, which is one order of magnitude 
smaller than the value reached in the numerical experiments of Lega et 
al. (2003). 

After r normalization steps, the Hamiltonian reads: 

= Z^(I, 4>) + 4i + +45 2 + •■■ (34) 

The normal form Z( r \l,<f>) contains all the Fourier terms up to \k\ < 
3r — 1 belonging to the resonant module Ai = {k : k\ — 2k<i = 
and &3 = 0}. The normal form Z^ r \l, <f>) is an integrable hamiltonian, 
the integrals, besides the energy, being /3 and If = 2ii + The level 
lines If = c are normal to the resonance line I\ — 2I2 = 0. Since 
we are interested in measuring the drift along the chaotic border of the 



+ 1 • 



(33) 
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Figure 4- (a) The size of the truncated remainder \\R^ r ' \\ a function of the 

truncation order r t > r after one and after two normalization steps (upper and lower 
points respectively) when e = 0.001. (b) The size of the remainder as a function of 
the order of the normalization for three different values of e. ||i? op t|| as a function 
of e is given by the size of the remainder at the minimum of each curve. 



resonance, we can estimate the values |ii| and II2I at the chaotic border 
either numerically, i.e., by measuring the half-width of the resonant 
zone in each panel of fig.2 of Lega et al. (2003), or by using Eq.(30) 
with k± = 1,&2 = —2. The numerical calculation was obtained by 
measuring the half- width AI of the interval along the line 2/i + 12 = c 
passing through (Iu,l2*) = (0.31,0.155) with endpoints at the two 
sides of the chaotic border marked yellow in Lega et al. (2003), and 
by dividing AI by e 1 / 2 according to the rescaling (18). As expected, 
the estimates for three different values of e turned to be approximately 
equal, yielding average values (|ii|, I/2I) = (0.0325, 0.065). These values 
were inserted in (34) in order to calculate the norm of the remainder 

R^ r \l,4>) = -R^+i + R r r +2 + •••• I* should be noticed that I3 does not 
appear in the remainder because it is a dummy action, i.e., it only 
appears linearly in the unperturbed Hamiltonian. 

According to standard theory, the remainder series R^{I,<p) should 
be convergent. This is shown in Fig.4a, in which the truncated sum 

ll^ (r) ll<^ £ Hitfll (35) 

j=r+l 

is plotted against the truncation order rt for the remainders of the first 
two normalization steps r = 1 and r = 2 for e = 0.001 (the norm 
1 1 • 1 1 is taken as the sum of moduli of all the trigonometric coefficients 
of the involved function). Clearly, the remainder shows the tendency 
to converge to its final size after summing just its first two or three 
successive terms. Furthermore, the size of the remainder, estimated as 
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||-R^||<15 exhibits the well known behavior expected for an asymptotic 
series (Fig. 4b). Namely, the size of the remainder initially decreases, as 
r increases, up to an optimal order r op t(e) corresponding to the minima 
of the curves of Fig. 4b for different values of e. Beyond the optimal 
order, however, the remainder becomes an increasing function of r and 
one has lim^oo ||i?( r )|| = oo, marking the eventual divergence of the 
normalization procedure. From fig. 5b we see clearly that the optimal 
order r op t increases as e decreases. 

Figures 5 and 6 show now the main result. The abscissa in Fig. 5 
shows the value of ||i? op i|| = ||-R^ r ° p *^^||<i5 for ten different values of 
e in the range 0.001 < e < 0.02. The ordinate shows the value of the 
diffusion coefficient D(e) for the same values of e as given by Lega et al. 
(2003). The straight line represents a power-law fitting of the relation 
between the optimal remainder and the diffusion coefficient. The best 
fit law is D = 1.02| \R op t 1 1 2 ' 97 . Thus, the scaling is essentially: 

D cc \\R op t\\ 3 ■ (36) 

On the other hand, Fig. 6 shows the scaling of ||i? pi|| with 1/e when 
e reaches values as small as e = 0.0001, i.e., one order of magnitude 
smaller than the last point of Fig. 7 of Lega et al. (2003) (the right ver- 
tical dashed line shows the last point where the numerical calculation 
of Lega et al. (2003) was stopped). The two solid lines passing through 
the data correspond to a power-law (lower curve) and an exponential 
law (upper curve) fitting the data. For small values of 1/e (i.e. for 
large values of e), the power-law fitting is better than the exponential 
fitting. The limit beyond which the exponential fitting is acceptable 
is around e = 0.01 (indicated by the left vertical dashed line of Fig.6 
at 1/e = 100). We may identify this limit as roughly corresponding to 
the threshold of the Nekhoroshev regime. On the other hand, in the 
interval 100 < 1/e < 1000, both the power law and the exponential law 
yield acceptable fittings. Nevertheless, beyond the value 1/e > 1000, 
the power-law clearly fails while the exponential fitting follows now 
narrowly the data. The numerical best fit exponential law yields 

\\Ropt\ \ oc exp(-e*/e a21 ) . 

We note that Froeschle et al. (2005), making numerical experiments 
of the diffusion of orbits in a mapping model resembling to the Hamilto- 
nian (22), found an exponential scaling of the diffusion coefficient with 
e in which the best-fit exponent is a = 0.28. However, in that paper 
too the exponential scaling was unraveled by considering values of e 
smaller by more than one orders of magnitude from a critical threshold 
value below which resonances do not significantly overlap. We conclude 
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Figure 5. The diffusion coefficient D reported in Lega et af. (2003) as a function of 
the size of the optimal remainder ||-R pt|| of the Birkhoff normal form calculated for 
different values of e. 



I 1 ' 




Figure 6. The size of the optimal remainder ||i? op t|| as a function of 1/e. The data 
clearly deviate from a power law beyond 1/e = 1000, and they are well fitted by an 
exponential law, as predicted by the Nekhoroshev theorem. The left vertical dashed 
line gives a lower threshold of 1/e below which the exponential law is no longer valid. 
The data are well fitted also by a power law in the range 10 < 1/e < 1000. The right 
vertical dashed line shows the last point of the numerical calculation of the diffusion 
coefficient by Lega et al. (2003) . 
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that the exponential law predicted by the Nekhoroshev theorem can 
be unambiguously seen in the real data by various methods when e 
becomes at least one order of magnitude smaller than the critical value 
characterizing the onset of the 'Nekhoroshev regime'. 

4. Conclusions and Discussion 

The present paper demonstrates the applicability of the analytical 
apparatus of the Nekhoroshev theory in recovering the exponential 
scaling of the coefficient of Arnold diffusion D = 0(exp(— 1/e)) along a 
particular resonance in the model Hamiltonian system of three degrees 
of freedom of Froeschle et al. (2000). In particular: 

1) The normalization algorithm for the Hamiltonian, based on an 
organization of the perturbation series in terms of the Fourier order of 
the terms (rather than on powers of e), is presented, and its benefits 
are analyzed in detail. 

2) The implementation of this algorithm on the computer allowed 
us to reach the optimal order of normalization, at which the size of the 
remainder ||i? op t|| of the normal form becomes minimal. 

3) The coefficient of Arnold diffusion D, as determined numerically 
by Lega et al. (2003), is found to scale with the size of the optimal 
remainder like D oc R^ pt . 

4) The size of the optimal remainder is found to scale exponen- 
tially with the inverse of the perturbation parameter, namely ||i? pt|| °c 
exp(— l/e a ), with a ~ 0.21. The exponential scaling clearly shows when 
one considers values of e as small as e = 10~ 4 , i.e., one order of 
magnitude smaller than in the numerical experiments of Lega et al. 
(2003). 

Estimates on the speed of Arnold diffusion based on numerical ex- 
periments have been given in the literature by a number of authors. 
We note in particular the early work of Kaneko and Konishi (1989), in 
which an exponential scaling law is found with a between 0.1 and 0.3, 
i.e., consistent with the present results. 

On the other hand, Wood et al. (1990) provided estimates of the 
diffusion coefficient in standard-like multidimensional symplectic maps 
on the basis of the Arnold-Melnikov method, i.e., the exponentially 
small splitting of separatrices due to the presence of higher order cou- 
pling terms in the resonant Hamiltonian. An exponential scaling of D 
with 1/e was also found in that case, favoring though a value of a 
rather close to a = 1/2. Thus, a detailed comparison of the results 
by the Nekhoroshev and Arnold-Melnikov theories is in order. To our 
knowledge, the only hint in that direction is a paper by Morbidelli and 
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Giorgilli (1997). Nevertheless, more work is necessary in order to clarify 
this connection by specific quantitative studies, as well as to compare 
the predictions of the two theories for Arnold diffusion in various types 
of multidimensional systems. 

Finally, the fact that the scaling of ||-R pt|| with e appears also as 
a power law in a transient interval of values of e, before the onset of 
the exponential regime, is consistent with some past theoretical work 
(Chirikov et al. 1985) calling such a transient regime 'modulational 
diffusion'. According to this theory, there is an intermediate interval 
of values of e within which many high order resonances, located in the 
chaotic border of the low-order resonance along which the diffusion 
takes place, still overlap. We conjecture that these are resonances of 
some order \k\ which must be above, but close to a truncation order K 
estimated as K ~ l/e a (according to the analysis of subsection 2.3). 

At any rate, the results of the present paper provide a clear relation 
between the local value of the diffusion coefficient D along a resonance, 
on the one hand, and the size of the optimal remainder ||i? op i|| of the 
normal form for the same resonance, on the other hand. Furthermore, 
for e sufficiently small ||i? pt|| is found to scale with e precisely as 
predicted by the Nekhoroshev theory. Thus, the final conclusion of our 
study that the analytical techniques of the Nekhoroshev theory can be 
used with much profit, in order to construct precise quantitative esti- 
mates of the speed of Arnold diffusion in Hamiltonian systems. 
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Appendix 

The following is a more technical analysis of point (b) of subsection 
2.4 referring to the organization of the perturbation series in powers of 
the small parameter e. 

Let 

r max be the maximum order in the book-keeping variable A = e, 
Kmax be the maximum order in Fourier space specifying a domain 
in which we are interested in computing the Hamiltonian, and K the 
order of truncation in Fourier space beyond which terms are not nor- 
malized (in general we have to set K max > K so that some remainder 
terms of the Hamiltonian be also stored). Figure 7a shows this do- 
main as a gray-shaded parallelogram in the space (r, \ k\). Now, it can 
be readily seen that many terms of the Hamiltonian which lie out- 
side the domain of interest at some particular order r interact with 
the generating function in such a way as to produce, at some sub- 
sequent order, terms which lie inside the domain of interest. Indeed, 
consider a term of the form e ri a( J) exp(ik ■ 4>) in the generating func- 
tion, with r\ < r max and \k\ < K < K rnax (inside the domain), 
and a term e r2 a'( J) ex.p(ik' • 0) in the Hamiltonian, with T2 < r max 
and | A; | > K max (outside the domain). Assume that r\ + ri < r max . 
Then, a Poisson bracket operation between the two terms yields a 
term e ri+r2 [a(J)k ■ Vja'(J) - a'(J)k' ■ Vja(J)] exp(i(k + k') ■ <f>). If 
\k + k'\ < K max (this is possible because the components of k and k' 
are added algebraically), the so produced term lies inside the domain 
of interest, despite the fact that the Hamiltonian term from which it 
originated lies outside the domain of interest. 

A careful examination of all the Lie operations taking place up 
to the order r max shows that in order to obtain complete knowledge 
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Figure 7. (a) The gray parallelogram shows the 'domain of interest' defined as the 
domain in Fourier space (orders of |fc|) of the terms that we wish to store in the 
transformed Hamiltonian up to an order r = r ma x in the book-keeping parameter 
A = e. The inclined line, given by \k\ = K max (r max — r + 1), represents the upper 
limit of a larger domain in Fourier space within which intermediate terms must be 
stored at an order r, yielding contributions in the gray-shaded domain at orders 
larger than r. (b) Under an appropriate book-keeping (see subsection (2.4)) the 
domain of interest can be made to coincide with the domain in which terms must 
be stored. 



of the transformed Hamiltonian in a domain {rmax ? ■^■max 

), one must 

store, for every order r < r max , all the Hamiltonian terms of Fourier 
order |fc| < K(r max — r) + K max . This extended domain, shown by the 
triangular domain of figure la, contains many more terms than those 
of the domain of interest (the ratio of the number of terms in the two 
domains is proportional to the ratio of the areas of the domains raised 
to a power ~ 2n, where n is the number of degrees of freedom). 

At this point we may argue that, for sufficiently large K max , the 
terms with \k\ > K max have very small size because of the analyticity 
condition (Eq.(3)), so that they can probably be ignored without great 
modification of the results. However, precisely this type of argument 
brings about the real source of the above problem, which is the choice of 
'book-keeping' A = e. The fact that the triangular domain of figure la 
contains many more terms than the gray parallelogram in the same 
figure simply depicts failure to recognize that for any order r < r max , 
many Fourier terms formally stored as of order A r are actually of much 
smaller size than e r . For example, while Hi is an overall 0(1) quantity, 
not all the Fourier terms in Hi have a similar size, because the size of an 
exp(ifc • (f>) term in Hi scales as (e~°")' fc L This initial scaling propagates 
also at subsequent orders of normalization. This fact suggests that the 
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terms should be 'book-kept' differently, i.e., on the basis of their Fourier 
order rather than their e-order. 

Following the algorithm exposed in subsection 2.4 (steps (1) to (3)), 
the domain of interest (figure lb) is made to coincide with the domain 
in which terms should be stored, i.e., terms within the gray-shaded 
area can only be produced by Poisson bracket operations between 
terms which are also within the gray shaded area. In particular, the 
meaning of step (2) is that in a Hamiltonian like (1) the parameter e is 
not so relevant in characterizing the smallness of the various terms in 
Hi (e can in fact be incorporated in the book-keeping scheme, as was 
done in section 3). A formal analogy can be made with a polynomial 
Hamiltonian around an elliptic equilibrium: 

H = j^ u Pl+A + e £ a mum2 _ mn q^qr-C n • (37) 

In the case of the Hamiltonian (37) it is customary to proceed by 
normalizing iteratively by the degrees of the monomial terms in the 
canonical variables rather than by powers of e. Such a choice under- 
lines that the real smallness in the polynomial case is the distance 
P = \fj2(<li + Pi) from the equilibrium point. In fact, it is often conve- 
nient to rescale the whole Hamiltonian (37) so that e disappears from 
the lowest order terms of the polynomial expansion in the r.h.s. of (37). 
The analogy with the generic case becomes clear by introducing action- 
angle variables % = \J1Ji sin (pi, pi = \/2Ji cos (pi, i = 1, . . . , n. We then 
readily find that a Fourier term of order p r ~ J r / 2 is always of Fourier 
order |fc| < r (and \k\ has the same parity as r). The relevant domain of 
interest is thus limited by the equation \k\ = r which yields a domain 
very similar to Fig. 7b (the only difference is in the slope of the limiting 
upper line, which is [l/er] in Fig.7b and 1 in the polynomial case). 
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